#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/Fiber.h"
#include "../DWMRICore/VectorField.h"

#include <vector>
#include <string>

enum SR_METHOD {
	SR_NONE,
	SR_FA,
	SR_TDI
};

char *pFiberConfigurePathname = NULL;
char *pDstPathname = NULL;
char *pFAPathname = NULL;
char *pDECPathname = NULL;

char **confPathname = NULL;
char **vertPathname = NULL;
int w, h, d, scale;
float stepSize, maxLength;
int res_scale = 4;

SR_METHOD method = SR_NONE;

float gaussianWeight[1100];

uint sub_steps = 5;


const float offsetX = 24.5f;
const float offsetY = 24.5f;
const float offsetZ = 24.5f;

void PrintUsage()
{
	printf("Usage: SuperResolutionVolume -fiber <fiber_configure_file> -dst <dst_file> -method <method> [-scale <scale> -FA <FA_file> -DEC <DEC_file>]\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-fiber") == 0) {
			pFiberConfigurePathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			pDstPathname = argv[++i];
		} else if (strcmp(argv[i], "-method") == 0) {
			++i;
			if (strcmp(argv[i], "FA") == 0) {
				method = SR_FA;
			} else if (strcmp(argv[i], "TDI") == 0) {
				method = SR_TDI;
			} else {
				method = SR_NONE;
			}
		} else if (strcmp(argv[i], "-scale") == 0) {
			res_scale = atoi(argv[++i]);
		} else if (strcmp(argv[i], "-FA") == 0) {
			pFAPathname = argv[++i];
		} else if (strcmp(argv[i], "-DEC") == 0) {
			pDECPathname = argv[++i];
		} else {
			return -1;
		}
	}

	if (pFiberConfigurePathname == NULL || pDstPathname == NULL || method == SR_NONE || res_scale == 0)
		return -2;

	if (method == SR_FA && pFAPathname == NULL)
		return -2;

	return 0;
}

float ComputeWeightedFA(CScalar *pFA, CPoint3F &v, std::vector<CPoint3F> &pointVec, CVector3F &dir)
{
	float fa = 0.0f;
	float weight = 0.0f;

	uint index = -1;
	float minDis = FLT_MAX;

	for (uint i = 0; i < pointVec.size(); ++i) {
		CVector3F d = pointVec.at(i) - v;
		float dis = d.Length();
		if (dis < minDis) {
			minDis = dis;
			index = i;
		}
	}

	uint prev_index = index < 1 ? 0 : (index-1);
	uint next_index = index >= (pointVec.size()-1) ? (pointVec.size()-1) : (index+1);
	dir = pointVec.at(next_index) - pointVec.at(prev_index);
	dir.Normalize();

	uint step = sub_steps * 2;
	for (uint i = 0; i < pointVec.size(); i+=step) {
		uint id;
		if (i > index)
			id = i - index;
		else
			id = index - i;

		if (id > 1000)
			continue;
		CPoint3F pt = pointVec.at(i);
		fa += gaussianWeight[id] * pFA->GetScalar(pt);
		weight += gaussianWeight[id];
	}

	return fa / weight;

	//return 1.0f;
}
//
//float ComputeWeightedFA(CScalar *pFA, CPoint3F &v, CPoint3F *pointBuffer, const uint size, CVector3F &dir)
//{
//	float fa = 0.0f;
//	float weight = 0.0f;
//
//	uint index = -1;
//	float minDis = FLT_MAX;
//
//	for (uint i = 0; i < size; ++i) {
//		CVector3F d = pointBuffer[i] - v;
//		float dis = d.Length();
//		if (dis < minDis) {
//			minDis = dis;
//			index = i;
//		}
//	}
//
//	uint prev_index = index < 1 ? 0 : (index-1);
//	uint next_index = index >= (size-1) ? (size-1) : (index+1);
//	dir = pointBuffer[next_index] - pointBuffer[prev_index];
//	dir.Normalize();
//
//	for (uint i = 0; i < size; i+=sub_steps) {
//		uint id;
//		if (i > index)
//			id = i - index;
//		else
//			id = index - i;
//
//		fa += gaussianWeight[id] * pFA->GetScalar(pointBuffer[i]);
//		weight += gaussianWeight[id];
//	}
//
//	return fa / weight;
//
//	//return 1.0f;
//}
//
void AddDirection(CVector3F &v, CVector3F dir)
{
	//if (v.Length() < 0.5f && v.Length() > 10e-7) {
	//	v = dir;
	//} else {
	//	if (InnerProduct(v, dir) > 0.0f) {
	//		v.m_x += dir.m_x; v.m_y += dir.m_y; v.m_z += dir.m_z;
	//	} else {
	//		v.m_x -= dir.m_x; v.m_y -= dir.m_y; v.m_z -= dir.m_z;
	//	}
	//}

	//if (InnerProduct(v, dir) > 0.0f) {
	//	v.m_x += dir.m_x; v.m_y += dir.m_y; v.m_z += dir.m_z;
	//} else {
	//	v.m_x -= dir.m_x; v.m_y -= dir.m_y; v.m_z -= dir.m_z;
	//}

	v.m_x += fabs(dir.m_x);
	v.m_y += fabs(dir.m_y);
	v.m_z += fabs(dir.m_z);
}
//
//void ComputeSR_FA()
//{
//	/* initialize the gaussian weight */
//	for (int i = 0; i < 1100; ++i) {
//		float dis = (float)i * stepSize / float(sub_steps);
//		gaussianWeight[i] = exp(-(dis * dis / (0.1f * maxLength * maxLength)));
//	}
//
//	CScalar *pFA = new CScalar();
//	pFA->OpenScalarFile(pFAPathname);
//
//	int sw = w * res_scale;
//	int sh = h * res_scale;
//	int sd = d * res_scale;
//	float *SRFA = new float[sw*sh*sd];
//	memset(SRFA, 0, sizeof(float)*sw*sh*sd);
//	uint *count = new uint[sw*sh*sd];
//	memset(count, 0, sizeof(uint)*sw*sh*sd);
//
//	CVector3F *SRDIR = new CVector3F[sw*sh*sd];
//	memset(SRDIR, 0, sizeof(CVector3F)*sw*sh*sd);
//
//	for (int i = 0; i < sh; ++i) {
//		printf("\ty = %d ... ", i);
//		CFiber *fibers;
//		int cw, ch;
//		ReadFibers(confPathname[i], vertPathname[i], &fibers, cw, ch);
//
//		CPoint3F *pointBuffer = new CPoint3F[1100];
//		uint point_count = 0;
//
//		for (int j = 0; j < cw * ch; ++j) {
//			if ((fibers[j].m_fCount + fibers[j].m_bCount) > 2) {
//				std::vector<int> indexVec;
//				memset(pointBuffer, 0, sizeof(CPoint3F)*1100);
//				point_count = 0;
//				for (uint k = fibers[j].m_fCount-1; k > 1 ; --k) {
//					CPoint3F s_pt = fibers[j].m_pF[k];
//					CPoint3F e_pt = fibers[j].m_pF[k-1];
//					for (uint m = 0; m < sub_steps; ++m) {
//						float ratio = (float)m / float(sub_steps);
//						CPoint3F pt = s_pt * (1.0f - ratio) + e_pt * ratio;
//						//pointVec.push_back(pt);
//						pointBuffer[point_count++] = pt;
//					}
//				}
//				for (uint k = 1; k < fibers[j].m_bCount; ++k) {
//					CPoint3F s_pt = fibers[j].m_pB[k-1];
//					CPoint3F e_pt = fibers[j].m_pB[k];
//					for (uint m = 0; m < sub_steps; ++m) {
//						float ratio = (float)m / float(sub_steps);
//						CPoint3F pt = s_pt * (1.0f - ratio) + e_pt * ratio;
//						//pointVec.push_back(pt);
//						pointBuffer[point_count++] = pt;
//					}
//				}
//				//pointVec.push_back(fibers[j].m_pB[fibers[j].m_bCount-1]);
//				pointBuffer[point_count++] = fibers[j].m_pB[fibers[j].m_bCount-1];
//
//				//for (uint k = 0; k < point_count; ++k) {
//				//	int cx = int(pointBuffer[k].m_x * (float)res_scale);
//				//	int cy = int(pointBuffer[k].m_y * (float)res_scale);
//				//	int cz = int(pointBuffer[k].m_z * (float)res_scale);
//
//				//	if (cx < 0 || cx >= sw || cy < 0 || cy >= sh || cz < 0 || cz >= sd)
//				//		continue;
//
//				//	int index = (cz * sh + cy) * sw + cx;
//				//	bool flag = false;
//				//	for (int m = 0; m < indexVec.size(); ++m) {
//				//		if (indexVec.at(m) == index)
//				//			flag = true;
//				//	}
//				//	if (false == flag) {
//				//		CPoint3F v;
//				//		v.m_x = ((float)cx/* + 0.5f*/) / (float)res_scale;
//				//		v.m_y = ((float)cy/* + 0.5f*/) / (float)res_scale;
//				//		v.m_z = ((float)cz/* + 0.5f*/) / (float)res_scale;
//				//		CVector3F dir;
//				//		//SRFA[index] += ComputeWeightedFA(pFA, v, pointVec, dir);
//				//		SRFA[index] += ComputeWeightedFA(pFA, v, pointBuffer, point_count, dir);
//				//		AddDirection(SRDIR[index], dir);
//				//		count[index]++;
//				//		indexVec.push_back(index);
//				//	}
//				//}
//
//				{
//					int cx = int(fibers[j].m_seed.pos.m_x * (float)res_scale);
//					int cy = int(fibers[j].m_seed.pos.m_y * (float)res_scale);
//					int cz = int(fibers[j].m_seed.pos.m_z * (float)res_scale);
//
//					if (cx < 0 || cx >= sw || cy < 0 || cy >= sh || cz < 0 || cz >= sd)
//						continue;
//
//					int index = (cz * sh + cy) * sw + cx;
//					bool flag = false;
//					for (int m = 0; m < indexVec.size(); ++m) {
//						if (indexVec.at(m) == index)
//							flag = true;
//					}
//					if (false == flag) {
//						CPoint3F v;
//						v.m_x = ((float)cx/* + 0.5f*/) / (float)res_scale;
//						v.m_y = ((float)cy/* + 0.5f*/) / (float)res_scale;
//						v.m_z = ((float)cz/* + 0.5f*/) / (float)res_scale;
//						CVector3F dir;
//						//SRFA[index] += ComputeWeightedFA(pFA, v, pointVec, dir);
//						SRFA[index] += ComputeWeightedFA(pFA, v, pointBuffer, point_count, dir);
//						AddDirection(SRDIR[index], dir);
//						count[index]++;
//						indexVec.push_back(index);
//					}
//				}
//
//			}
//		}
//		SafeDeleteArray(fibers);
//		SafeDeleteArray(pointBuffer);
//
//		printf("done.\n");
//	}
//
//	for (int i = 0; i < sw*sh*sd; ++i) {
//		if (count[i] != 0) {
//			SRFA[i] = SRFA[i] / (float)count[i];
//			SRDIR[i] = SRDIR[i] / (float)count[i];
//			SRDIR[i] = SRDIR[i] * SRFA[i];
//		}
//	}
//
//	CScalar *pSR = new CScalar();
//	pSR->CreateScalar(sw, sh, sd, SRFA);
//	pSR->SaveScalarFile(pDstPathname);
//
//	CVectorField *pSRDIR = new CVectorField();
//	pSRDIR->CreateVectorField(sw, sh, sd);
//	for (int i = 0; i < sw*sh*sd; ++i) {
//		pSRDIR->SetVector(i, SRDIR[i]);
//	}
//	pSRDIR->SaveVectorFieldFile(pDECPathname);
//
//	SafeDeleteArray(SRFA);
//	SafeDeleteArray(SRDIR);
//	
//	SafeDelete(pFA);
//	SafeDelete(pSR);
//	SafeDelete(pSRDIR);
//}
//
//void MappingFiberTDI(const CFiber *fiber, float *TDI, int *mask, const int fiberID)
//{
//	int sw = w * res_scale;
//	int sh = h * res_scale;
//	int sd = d * res_scale;
//
//	//char *temp = new char[sw*sh*sd];
//	//memset(temp, 0, sizeof(char)*sw*sh*sd);
//
//	for (uint i = 1; i < fiber->m_fCount; ++i) {
//		CPoint3F p1 = fiber->m_pF[i];
//		CPoint3F p0 = fiber->m_pF[i-1];
//		CVector3F v = (p1 - p0) / 10.0f;
//		
//		CPoint3F q = p0;
//		for (int k = 0; k < 11; ++k) {
//			q = q + v;
//			
//			int cx = int(floor(q.m_x * (float)res_scale));
//			int cy = int(floor(q.m_y * (float)res_scale));
//			int cz = int(floor(q.m_z * (float)res_scale));
//
//			if (cx >= 0 && cx < sw && cy >= 0 && cy < sh && cz >= 0 && cz < sd) {
//				if (mask[(cz*sh+cy)*sw+cx] < fiberID)
//					TDI[(cz*sh+cy)*sw+cx] += 1.0f;
//				mask[(cz*sh+cy)*sw+cx] = fiberID;
//			}
//		}
//	}
//
//	for (uint i = 1; i < fiber->m_bCount; ++i) {
//		CPoint3F p1 = fiber->m_pB[i];
//		CPoint3F p0 = fiber->m_pB[i-1];
//		CVector3F v = (p1 - p0) / 10.0f;
//		
//		CPoint3F q = p0;
//		for (int k = 0; k < 11; ++k) {
//			q = q + v;
//			
//			int cx = int(floor(q.m_x * (float)res_scale));
//			int cy = int(floor(q.m_y * (float)res_scale));
//			int cz = int(floor(q.m_z * (float)res_scale));
//
//			if (cx >= 0 && cx < sw && cy >= 0 && cy < sh && cz >= 0 && cz < sd) {
//				if (mask[(cz*sh+cy)*sw+cx] < fiberID)
//					TDI[(cz*sh+cy)*sw+cx] += 1.0f;
//				mask[(cz*sh+cy)*sw+cx] = fiberID;
//			}
//		}
//	}
//
//	/*for (int i = 0; i < sw*sh*sd; ++i) {
//		if (mask[i] == fiberID) {
//			TDI[i] += 1.0f;
//		}
//	}*/
//
//	//SafeDelete(temp);
//}
//
//void ComputeSR_TDI()
//{
//	int sw = w * res_scale;
//	int sh = h * res_scale;
//	int sd = d * res_scale;
//	float *TDI = new float[sw*sh*sd];
//	memset(TDI, 0, sizeof(float)*sw*sh*sd);
//
//	for (int z = 0; z < sd; ++z) {
//		CFiber *fibers;
//		int cw, ch;
//		ReadFibers(confPathname[z], vertPathname[z], &fibers, cw, ch);
//
//		int *mask = new int[sw*sh*sd];
//		memset(mask, 0, sizeof(int)*sw*sh*sd);
//
//		for (int j = 0; j < cw * ch; ++j) {
//			if ((fibers[j].m_fCount + fibers[j].m_bCount) > 2) {
//				// mapping fiber to the super resolution volume
//				MappingFiberTDI(&(fibers[j]), TDI, mask, j);
//			}
//			float ff = float(j+1) / float(cw*ch) * 100.0f;
//			printf("\r%.2f%% ", ff);
//		}
//		delete[] fibers;
//		SafeDelete(mask);
//
//		printf("fibers on z = %d processed.\n", z);
//	}
//
//	CScalar *pSR = new CScalar();
//	pSR->CreateScalar(w*res_scale, h*res_scale, d*res_scale, TDI);
//	pSR->SaveScalarFile(pDstPathname);
//	delete pSR;
//}
//
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	/* read the configure file */
//	ReadConfigureFile(pFiberConfigurePathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);
//	//w = w / scale;
//	//h = h / scale;
//	//d = d / scale;
//
//	int temp = d;
//	d = h;
//	h = temp;
//
//	w = w / scale;
//	h = h / scale;
//	d = d / scale;
//
//	if (method == SR_FA) {
//		/* compute super resolution FA volume */
//		ComputeSR_FA();
//	} else if (method == SR_TDI) {
//		/* compute super resolution TDI volume */
//		ComputeSR_TDI();
//	}
//
//	return 0;
//}


int FindOffset(FILE *fp)
{
	char buffer[256];
	do {
		fscanf(fp, "%s", buffer);
	} while (strcmp(buffer, "file:") != 0);

	fscanf(fp, "%s", buffer);
	
	int offset;
	fscanf(fp, "%d", &offset);

	return offset;
}

int FindFiberCount(FILE *fp)
{
	char buffer[256];
	do {
		fscanf(fp, "%s", buffer);
	} while (strcmp(buffer, "count:") != 0);

	int count = 0;
	fscanf(fp, "%d", &count);

	return count;
}

void ComputeSR_TDI(const int size)
{
	int sw = w * res_scale;
	int sh = h * res_scale;
	int sd = d * res_scale;
	float *SRTDI = new float[sw*sh*sd];
	memset(SRTDI, 0, sizeof(float)*sw*sh*sd);

	for (int i = 0; i < size; ++i) {
		FILE *fp = fopen(vertPathname[i], "rb");

		/* find the offset of the binary data */
		int offset = FindOffset(fp);
		/* find the fiber count */
		fseek(fp, 0, SEEK_SET);
		int fiberCount = FindFiberCount(fp);
		
		if (fseek(fp, offset, SEEK_SET) == 0) {
			for (int k = 0; k < fiberCount; ++k) {
				/* read the fiber */
				std::vector<CPoint3F> fiber;
				float p[3];
				while(true) {
					fread(p, sizeof(float), 3, fp);
					if (p[0] != p[0])
						break;
					//CPoint3F pt((p[0]+32.5f)/0.507812f, (p[1]+32.5f)/0.507812f, (p[2]+17.5f)/0.5f);
					CPoint3F pt((p[0]+offsetX), (p[1]+offsetY), (p[2]+offsetZ));
					fiber.push_back(pt);
				}

				/* resample the fiber */
				std::vector<CPoint3F> pts;
				for (uint n = 1; n < fiber.size(); ++n) {
					CPoint3F s_pt = fiber.at(n-1);
					CPoint3F e_pt = fiber.at(n);
					for (uint m = 0; m < sub_steps; ++m) {
						float ratio = (float)m / float(sub_steps);
						CPoint3F pt = s_pt * (1.0f - ratio) + e_pt * ratio;
						pts.push_back(pt);
					}
				}
				pts.push_back(fiber.back());

				std::vector<int> indexVec;
				/* mapping to TDI */
				for (uint n = 0; n < pts.size(); ++n) {
					int cx = int(pts.at(n).m_x * (float)res_scale);
					int cy = int(pts.at(n).m_y * (float)res_scale);
					int cz = int(pts.at(n).m_z * (float)res_scale);

					if (cx < 0 || cx >= sw || cy < 0 || cy >= sh || cz < 0 || cz >= sd)
						continue;

					int index = (cz * sh + cy) * sw + cx;
					bool flag = false;
					for (int m = 0; m < indexVec.size(); ++m) {
						if (indexVec.at(m) == index) {
							flag = true;
							break;
						}
					}
					if (false == flag) {
						SRTDI[index] += 1.0f;;
						indexVec.push_back(index);
					}
				}
			}
		}
		
		fclose(fp);
	}

	CScalar *pSR = new CScalar();
	pSR->CreateScalar(sw, sh, sd, SRTDI);
	pSR->SaveScalarFile(pDstPathname);


	SafeDelete(SRTDI);
}

void ComputeSR_FA(const int size)
{
	/* initialize the gaussian weight */
	for (int i = 0; i < 1100; ++i) {
		float dis = (float)i * stepSize / float(sub_steps);
		gaussianWeight[i] = exp(-(dis * dis / (0.1f * maxLength * maxLength)));
	}

	int sw = w * res_scale;
	int sh = h * res_scale;
	int sd = d * res_scale;
	float *SRFA = new float[sw*sh*sd];
	memset(SRFA, 0, sizeof(float)*sw*sh*sd);
	uint *count = new uint[sw*sh*sd];
	memset(count, 0, sizeof(uint)*sw*sh*sd);

	//CVector3F *SRDIR = new CVector3F[sw*sh*sd];
	//memset(SRDIR, 0, sizeof(CVector3F)*sw*sh*sd);

	CScalar *pFA = new CScalar();
	pFA->OpenScalarFile(pFAPathname);


	for (int i = 0; i < size; ++i) {
		FILE *fp = fopen(vertPathname[i], "rb");

		/* find the offset of the binary data */
		int offset = FindOffset(fp);
		/* find the fiber count */
		fseek(fp, 0, SEEK_SET);
		int fiberCount = FindFiberCount(fp);

		if (fseek(fp, offset, SEEK_SET) == 0) {
			for (int k = 0; k < fiberCount; ++k) {
				/* read the fiber */
				std::vector<CPoint3F> fiber;
				float p[3];
				while(true) {
					fread(p, sizeof(float), 3, fp);
					if (p[0] != p[0])
						break;
					CPoint3F pt((p[0]+32.5f)/0.507812f, (p[1]+32.5f)/0.507812f, (p[2]+17.5f)/0.5f);
					//CPoint3F pt((p[0]+offsetX), (p[1]+offsetY), (p[2]+offsetZ));
					fiber.push_back(pt);
				}

				/* resample the fiber */
				std::vector<CPoint3F> pts;
				for (uint n = 1; n < fiber.size(); ++n) {
					CPoint3F s_pt = fiber.at(n-1);
					CPoint3F e_pt = fiber.at(n);
					for (uint m = 0; m < sub_steps; ++m) {
						float ratio = (float)m / float(sub_steps);
						CPoint3F pt = s_pt * (1.0f - ratio) + e_pt * ratio;
						pts.push_back(pt);
					}
				}
				pts.push_back(fiber.back());

				/* SRFA */
				std::vector<int> indexVec;
				for (uint k = 0; k < pts.size(); ++k) {
					int cx = int(pts.at(k).m_x * (float)res_scale);
					int cy = int(pts.at(k).m_y * (float)res_scale);
					int cz = int(pts.at(k).m_z * (float)res_scale);

					if (cx < 0 || cx >= sw || cy < 0 || cy >= sh || cz < 0 || cz >= sd)
						continue;

					int index = (cz * sh + cy) * sw + cx;
					bool flag = false;
					for (int m = 0; m < indexVec.size(); ++m) {
						if (indexVec.at(m) == index)
							flag = true;
					}
					if (false == flag) {
						CPoint3F v;
						v.m_x = ((float)cx/* + 0.5f*/) / (float)res_scale;
						v.m_y = ((float)cy/* + 0.5f*/) / (float)res_scale;
						v.m_z = ((float)cz/* + 0.5f*/) / (float)res_scale;
						CVector3F dir;
						SRFA[index] += ComputeWeightedFA(pFA, v, pts, dir);
						//SRFA[index] += ComputeWeightedFA(pFA, v, pointBuffer, point_count, dir);
						//AddDirection(SRDIR[index], dir);
						count[index]++;
						indexVec.push_back(index);
					}
				}
			}
		}

		fclose(fp);
	}

	//for (int i = 0; i < sw*sh*sd; ++i) {
	//	if (count[i] != 0) {
	//		SRFA[i] = SRFA[i] / (float)count[i];
	//	}
	//}

	std::string SRFAPathname;
	SRFAPathname = pDstPathname;
	SRFAPathname += "_SRFA.nrrd";

	CScalar *pSR = new CScalar();
	pSR->CreateScalar(sw, sh, sd, SRFA);
	pSR->SaveScalarFile(SRFAPathname.c_str());

	std::string CountPathname;
	CountPathname = pDstPathname;
	CountPathname += "_COUNT.nrrd";

	CScalar *pCount = new CScalar();
	pCount->CreateScalar(sw, sh, sd);
	for (int i = 0; i < sw*sh*sd; ++i) {
		if (count[i] == 0) {
			pCount->SetScalar(i, 1.0f);
		} else {
			pCount->SetScalar(i, (float)(count[i]));
		}
	}
	pCount->SaveScalarFile(CountPathname.c_str());

	
	//CVectorField *pSRDIR = new CVectorField();
	//pSRDIR->CreateVectorField(sw, sh, sd);
	//for (int i = 0; i < sw*sh*sd; ++i) {
	//	pSRDIR->SetVector(i, SRDIR[i]);
	//}
	//pSRDIR->SaveVectorFieldFile(pDECPathname);


	SafeDeleteArray(count);
	SafeDeleteArray(SRFA);
	//SafeDeleteArray(SRDIR);
	

	SafeDelete(pFA);
	SafeDelete(pSR);
	//SafeDelete(pCount);
	//SafeDelete(pSRDIR);
}


int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	/* read the configure file */
	FILE *fp = fopen(pFiberConfigurePathname, "r");
	fscanf(fp, "%d %d %d", &w, &h, &d);
	fscanf(fp, "%f", &stepSize);
	int count = 0;
	fscanf(fp, "%d", &count);
	vertPathname = new char*[count];
	for (int i = 0; i < count; ++i) {
		vertPathname[i] = new char[512];
		fscanf(fp, "%s", vertPathname[i]);
	}

	fclose(fp);

	maxLength = 10.0f;

	if (method == SR_TDI) {
		/* compute super resolution TDI volume */
		ComputeSR_TDI(count);
	} else if (method == SR_FA) {
		ComputeSR_FA(count);
	}

	return 0;
}